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ABSTRACT 

We estimate the two- and three-dimensional power spectra, P2(K) and Ps(k), of the 
galaxy distribution by applying a maximum likelihood estimator to pixel maps of the 
APM Galaxy Survey. The analysis provides optimal estimates of the power spectra 
and of their covariance matrices if the fluctuations are assumed to be Gaussian. Our 
estimates of P2{K) and Pa(k) are in good agreement with previous work but we 
find that the errors at low wavenumbers have been underestimated in some earlier 
studies. If the galaxy power spectrum is assumed to have the same shape as the 
mass power spectrum, then the APM maximum likelihood Pa(k) estimates at k < 
0.19/iMpc~ constrain the amplitude and shape parameter of a scale-invariant CDM 
model to lie within the 2a ranges 0.78 < (a s ) g < 1.18 and 0.05 < T < 0.38. Using the 
Galactic extinction estimates of Schlegel, Finkbeiner and Davis, we show that Galactic 
obscuration has a negligible effect on galaxy clustering over most of the area of the 
APM Galaxy Survey. 



1 INTRODUCTION 

In this paper we analyse the three-dimensional power spec- 
trum of galaxy clustering using the APM Galaxy Survey 
(Maddox et al. 1990a, b, c). The APM Galaxy Survey 
is a two-dimensional catalogue of galaxies complete to a 
magnitude limit of bj = 20.5 and covering an area of ap- 
proximately 12 percent of the sky. The survey has been 
used to estimate the angular two-point correlation function 
w(6) and the angular power spectrum P2(K), which are re- 
lated to their three-dimensional analogues £(r) and P3(k) via 
simple integral equations (Limber 1953; Groth and Peebles 
1977; Baugh and Efstathiou 1994). Recovering the three- 
dimensional power spectrum from angular statistics there- 
fore requires stable numerical techniques for inverting these 
integral equations. 

Baugh and Efstathiou (1993) described a technique for 
recovering the three-dimensional power spectrum from mea- 
surements of the angular correlation function. The three- 
dimensional power spectrum was parameterized by a set of 
amplitudes P3 1 (or 'bandpowers') over bands of wavenum- 
bers centred at wavenumber ki. The integral equation relat- 
ing w(8) to Pa(k) was solved using Lucy's (1974) iterative 
deconvolution technique. A similar technique was applied 
by Baugh and Efstathiou (1994) to recover Pa(fc) from es- 
timates of the two-dimensional power spectrum P2(K) and 
by Baugh (1996) to recover £(r) from w(0). These inves- 
tigations show that Lucy's algorithm can provide a stable 
inversion. However, it is difficult to derive a reliable covari- 
ance matrix for the recovered estimates of P|. Baugh and 



Efstathiou (1993, 1994) derived estimates of the errors by 
computing the scatter in the P3 derived from four nearly 
equal areas of the APM survey. However, since the number 
of zones is small, these error estimates are crude and cannot 
be used to fit theoretical models with any confidence. 

Recently, Dodelson and Gaztanaga (2000) have de- 
scribed a method of inverting w(8) to recover P3 that em- 
ploys a Bayesian prior to contrain the smoothness of the 
inversion. This method can return a covariance matrix for 
P3, but requires an estimate of the covariance matrix of the 
input estimates of w(6i) and a model for the Bayesian prior. 
Eisenstein and Zaldarriaga (1999) present another inversion 
technique using singular value decomposition (see e.g. Press 
et al. 1992). Their method also recovers the covariance ma- 
trix for P3 but requires estimates of w(6) and its covariance 
matrix as inputs. 

The purpose of this paper is two- fold. Firstly, to assess 
the effects of Galactic extinction on large scale clustering 
in the APM Survey using the extinction model of Schlegel, 
Finkbeiner and Davis (1998, hereafter SFD) based on the 
COBE/DIRBE and IRAS maps. Secondly, to apply to the 
APM Survey modern maximum likelihood (ML) techniques 
similar to those used to estimate the power spectrum of the 
cosmic microwave background (CMB) anisotropics (Bond, 
Jaffe and Knox 1998; de Bernardis et al. 2000; Hannay et 
al., 2000). With the increase in computer power over the ten 
years since the APM survey was completed, it is now feasi- 
ble to perform a direct ML estimate of the angular power 
spectrum over wavenumbers extending into the non-linear 
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regime. This provides an optimal estimate (under certain 
assumptions) of the power spectrum and its covariance ma- 
trix in a conceptually straightfoward way, avoiding the need 
for estimators of P2(K) or w(8) that require a model of the 
true power spectrum. (See e.g. Hamilton, 1997a, b; Tegmark 
1997, Kerscher et al. 2000, and references therein for a dis- 
cussion of estimators of P2(K) and w(8)). An additional 
advantage of ML methods is that it is as easy to compute 
bandpower estimates of the three-dimensional power spec- 
trum P 3 (and its covariance matrix) as it is to estimate the 
two-dimensional power spectrum. The inversion from two 
to three dimensions can therefore be done with the same 
computer code and without the need for any assumptions 
other than that the underlying fluctuations obey Gaussian 
statistics. 

The outline of this paper is as follows. Section 2 de- 
scribes the method and applies it to Gaussian realizations 
of the APM Survey. In Section 3, we use the SFD dust maps 
to show how the two-dimensional power spectrum is affected 
by Galactic extinction. A model for the mean distribution 
of galaxies with redshift dN(z)/dz is constructed using data 
from the 2dF Galaxy Redshift Survey and this is used to 
compute the two- and three-dimensional power spectra by 
ML. Constraints on theoretical models are discussed in Sec- 
tion 4 and our conclusions are summarized in Section 5. 



2 METHOD 

2.1 Relations between power spectra and 
correlation functions 

In this Section we follow the notation of Baugh and Efs- 
tathiou (1993, 1994, hereafter refered to as BE93 and BE94). 
The angular correlation function w(8) is related to the spa- 
tial correlation function £(r, t) via the relativistic form of 
Limber's equation 



w(6) 



2 JT JT x 4 F- 2 a 6 p 2 (x)t{r,t)dxdu 
[J -00 x 2 F~ 1 a 3 p(x)dxj 2 



(1) 



Peebles (1980, §50.16). In this equation, p(x) is the selec- 
tion function of the survey (the probability that a galaxy 
at coordinate distance x is detected in the survey), a is the 
cosmological scale factor, and the metric is 



ds 



c 2 dt 2 — a 2 [dx 2 /F 2 (x) + x 2 dd 2 + x 2 sin 2 8 > 



(2) 



Equation (|l|) assumes that the clustering of galaxies is in- 
dependent of luminosity. However, this is quite a weak as- 
sumption for a magnitude limited optical survey since most 
of the galaxies have luminosities in a narrow range around 
the characteristic luminosity C of the Schechter (1976) lu- 
minosity function. The physical separation between galaxy 
pairs separated by an angle 8 on the sky is 



a \ u/F*{x) + x*6 



2.2-1 1/2 



(3) 



where we have assumed that the angle 8 is small. In the rest 
of this paper we adopt a spatially flat cosmological model 
with matter density parameter Q m =0.3 and a cosmological 
constant contributing Qa ~ 0.7. 

The spatial correlation function £(r,t) is related to the 
three-dimensional power spectrum P 3 (k,i) by 



„ ,, , sin(fcr/a) , 2 ,, 

p 3 (k,t) > k 2 dk, 

(kr/a) 



(4) 



and following BE93 we will assume that Pz[k, t) is a sepa- 
rable function of comoving wavenumber k and redshift z. 



P 3 (k,t) 
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The two-dimensional power spectrum Ps(K) is related 
to the angular correlation function by 

{' 30 

P 2 (K) = 2it w(8)J (K8)8d8. (6) 
Jo 

From equations (Q), (^)-@, the two-dimensional power 
spectrum is related to the three dimensional power spec- 
trum by the integral equation 



KP 2 (K)= / g(K/k)P 3 (k)dk, 
Jo 



where the kernel g(K/k) is 



g(K/k) 



M 2 Q. 2 



(f)' 



( dz\ 2 

\dx) (1 



F(x) 
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(7a) 



(7b) 



x=K/k 



(see BE94) and we have written the selection function p(x) 
in terms of the redshift distribution dN/ dz of the sample 



p(x) 



1 F(x) dN dz 



(8) 



M 2 Q 2 x 2 a? dz dx ' 

where J\f is the mean surface density of galaxies and f2 s is 
the solid angle of the survey. If we know the redshift distri- 
bution of a two-dimensional survey, the three-dimensional 
power spectrum can be recovered from estimates of the 
two-dimensional power spectrum by inverting equation ( [7a] ) 
using, for example, Lucy's (1974) method as described by 
BE94. However, in the next section we show that the inver- 
sion can be done by using a ML estimator. The ML method 
actually solves two problems simultaneously, solving the in- 
version probem and providing an optimal estimator of the 
power spectra Pz(K) and P 3 (k). 



2.2 Maximum likelihood estimator 

Assume that the galaxy catalogue is pixelized into a map of 
N identical pixels with galaxy count m in the i'th pixel. We 
define the data vector A as 



rii — (n) 



(9) 



where (n) is the mean galaxy count per pixel. 

If we assume that the A, constitute a Gaussian random 
field, the likelihood function is 



£ = 



— ±— 7 ^exp(-iA T C- 1 A 



{2n) N / 2 \C 
where C is the covariance matrix 
dj = (AiAA . 
From the definition of A;, 
1 



(10a) 



(10b) 



(n) 



Ca = -r-rSij + w(8ij), (11a) 
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where for square pixels of width 6 C 

®(*«) = 2^ XT Jo P ^ K ) K ° OS ( K9 ^ 

xW*(K6 c ,<f>)d<j>dK, 

and 



W C (K0 C ,4>) 



K8 C 



K6 C 



\ 2 



( 



(lib) 



(11c) 



For angul ar se parations much greater than the pixel size, 
equation (lib) simplifies to. 



P»{K)KJa{KOa)dK, 



> i 



(12) 



Equations (lib) and ( |l^ ) have been derived in the small 
angle limit 9j <C 1, which is a good approximation for the 
APM Galaxy Survey. This assumption is easily dropped, 
however, in which case equation (fl2|) reads 



Wis) = h EtW + l)ft(/)P£(coefl«). 



(13) 



In analogy with analyses of cosmic microwave background 
anisotropies, the angular wavenumber K is equivalent to the 
multipole moment I and the angular power spectrum P2(K) 
is equivalent to Ce (see e.g. Bond and Efstathiou 1987). 

Follo wing Bond, Jaffe and Knox (1998), the likelihood 
function (10a) can be maximized iteratively with respect to 
a set of parameters a p . Starting from an initial guess for the 
a p , the changes in the parameters 5a p at each iteration are 
calculated from 



5a„ = 



5S 



F~]Tt 

vv 



(AA ~C)[C 



dC 
da„i 



C~ 



where F„„i is the Fisher matrix 



d 2 \nC 

da p da p i 



-Tr 



c -idC_ c -i dC 



do 



(14) 



(15) 



The parameters a p can be chosen to be bandpower es- 
timates of the two-dimensional power spectrum P%(K) or 
of the three-dimensional power-spectrum Ps(k). Fo r th ese 
cases, the angular correlation function in equation (11a) is 
computed from the sum 



p 



where 



i 

2tt- 



1 
2tt- 



Jk P p +1 Jo Kcos (KOijCost) 

xW?(K6 c ,(j>) d(j>dK, 
' JT / fc fc ; +1 Jo 3iK/k) cos (KOacoH 

xwf(K6 c ,<j)) d<j>dkdK, 



(16) 

for P 2 (K), 
for P 3 (k). 



These integrals depend only on the binning of the parame- 
ters and on the pixel scale, so they can be computed once 
and stored. The computing time required to find the ML is 
dominated by the computation of the inverse matrix C _1 
and the multiplication of N x N matrices (both of which 
scale as TV 3 ). Our implementation on an SGI Origin 200 
workstation takes a few hours to converge to a solution for 
N « 4000. 



2.3 Tests of the Method 

We have tested the algorithm on simulated two-dimensional 
Gaussian random fields. We assume that the three- 
dimensional power spectrum is that of a linear adiabatic 
scale-invariant CDM model with a shape parameter of T = 
0.2 in the parameterization of Efstathiou, Bond and White 
(1992). The two-dimensional power spectrum was computed 
from equation ( |7a| ) using a model for the redshift distribu- 
tion of the APM Survey limited at b.j — 19.5 (see Section 
3.2 below). We adopt an evolution parameter of a — and 
normalize the spectra so that the rms fluctuation ampli- 
tude of the galaxy distribution averaged in spheres of radius 
8ft _1 Mpc spheres, (a 8 ) S) is unity. We used an 1024 2 FFT 
to generate a periodic Gaussian density field from the two- 
dimensional power spectrum in a 400° x 400° square from 
which we selected a 100° x 100° patch regridded into 32 x 32 
pixels for input into the ML code. The pixel size of the in- 
put catalogues is therefore 3.12°, but they include small scale 
power because they were generated on a grid of much higher 
resolution. 

The ML reconstructions averaged over 40 simulations 
are shown in Figure 1. Convergence to the ML solution for 
both the two- and three-dimensional power spectra is usually 
achieved within 5-10 iterations. The error bars shown on the 
points are computed from the inverse of the Fisher matrix, 
erf = F^ 1 , and are in excellent agreement with the scatter 
between simulations. 

There are a few subtle points about the analysis worth 
some discussion: 

[1] The sums over the bandpower parameters a p in equa- 
tion (]l^) are performed over a finite range of wavenumber 
K m in < K < K m ax (or k m i„ < fc < k m ax, depending on 
whether we are estimating the two- or three-dimensional 
power spectra). Ignoring power from wavenumbers outside 
these ranges leads to small biases in the ML solutions. 
In the examples shown in Figure 1, we have explicitly in- 
cluded integrals over the power spectra at K < K m in and 
K > Kmax assuming the input target power spectrum which 
is, of course, known. This removes any biases at large and 
small wavenumbers as shown in Figure 1. In application to 
real data, the power spectrum is unknown. In this case, one 
can simply increase the number of parameters extending the 
range of Kmin and K max and marginalize over a small num- 
ber (one should suffice) of parameters at either end of the 
wavenumber range. The remaining parameters will then be 
free of any biases. 

[2] The pixel scale of the maps used to generate Figure 1 
corresponds to a wavenumber K = 2n /9 C ~ 125. Neverthe- 
less, by correctly including t he wi ndow function of the pixels 
in the integral of equation ( |llb| ), the power spectrum can 
be recovered free of bias on sub-pixel scales, but obviously 
the errors become large as the estimates are extended below 
the pixel scale. In our application to the APM Survey, the 
limit on the pixel size is set by size of the data vector that 
can be analyzed in a reasonable amount of computer time. 
We find that it is possible to analyze maps with pixels of 
size 8 C = 0.89° (N w 4000 pixels) easily using workstations. 
It would be possible to increase the number of pixels by 
using supercomputers and by using Monte-Carlo methods 
as described by Oh, Spergel and Hinshaw (1999). However, 
in the ML analysis it is assumed that the underlying den- 
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Figure 1. ML bandpower estimates of the power spectra determined from simulated Gaussian random fields generated with a CDM 
power spectrum with shape parameter T = 0.2. The simulations mimic a square patch of the sky of side 100° with the selection function 
of the APM Galaxy Survey limited to bj = 19.5. The solid lines show the input two- and three-dimensional power spectra. The points 
show the mean bandpower estimates computed from 40 maps each containing 1024 pixels. The vertical error bars show lcr error estimates 
for a single simulation derived from the Fisher matrix. The horizontal error bars show the widths of each band. 



sity fluctuations are Gaussian, whereas the galaxy distribu- 
tion is observed to be strongly non-Gaussian on small scales 
where the distribution is also non-linear. At magnitude lim- 
its of bj « 19.5 — 20.0, the angular scales of significant 
non-Gaussianity and non-linearity in the APM survey are 
at K ^ 200. The ML estimator is therefore not guaranteed 
to be optimal or even unbiased at wavenumbers higher than 
K 200. This differs from the case of applying ML to the 
CMB anisotropies, where the assumption of Gaussian fluc- 
tuations is physically reasonable for primary anisotropies on 
all angular scales. 

[31 The numerical inversion of an integral equation such as 
( |7a[ ) is unstable; the inverted Ps(k) can show wild fluc- 
tuations as the number of bandpowers is increased (see 
e.g. BE93, BE94; Dodelson and Gaztanaga 2000). The ML 
method described here imposes no constraints on the band- 
power estimates a p and so there is no guarantee that the 
recovered power spectra will be smooth. As the number of 
bandpowers is increased, the ML solutions (particularly for 
Ps(k)) will begin to show oscillations. However, if we fit a 
theoretical model characterized by a few parameters to the 
data using the full covariance matrix of the estimates (see 
Section 4), then the best fitting parameters will be insen- 
sitive to the number of bandpowers and to oscillations in 
P 3 (k). 

[4] In analyzing the simulations, the mean galaxy count per 
pixel was estimated from each map by computing 



This is not strictly correct, since (n) is the mean galaxy 
count averaged over an ensemble of catalogues not the mean 
pixel count of a single map. This can introduce a bias that 
is related to the 'integral constraint' bias in estimates of 
w(ff) (e.g. Groth and Peebles 1977) and the power spectrum 
(e.g. Tadros and Efstathiou 1996). More correctly, the mean 
galaxy count should be treated as a para meter in the likeli- 
hood analysis. Maximising the likelihood (10a) with respect 
to (n) gives 



(n) 



E !7 ^iCr 1 



(18) 



(17) 



and so depends on the ML solution for the power spectrum. 
In practice the APM Galaxy Survey covers a large enough 
area that any bias introduced in using equation (hjj) is neg- 
ligible. 



3 APPLICATION TO THE APM GALAXY 
SURVEY 

In Section 3.1, we discuss the effects of Galactic extinction 
in the APM Galaxy Survey using the SFD dust map. This 
allows us to delineate an area of the APM Survey in which 
extinction has a negligible effect on the power spectrum. In 
Section 3.2, we use results from a small subset of the 2dF 
Galaxy Redshift Survey (see e.g. Colless, 1999) to derive 
a model for the redshift distribution of the APM Survey, 
improving on the model used by BE93, BE94. Results from 
the ML method are presented in Section 3.3. 
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Figure 2. Equal area projections of Galactic extinction and various galaxy maps. Clockwise from the top left: Extinction map in the 
southern Galactic cap computed from SFD; the extended APM Galaxy Survey limited at bj = 20; the APM Survey in the region with 
dust extinction of less than 0.2 magnitudes; the APM Survey in the region with dust extinction of less than 0.1 magnitudes. 



3.1 Input APM Galaxy Catalogue 

The APM Galaxy Survey is described in detail in a series of 
papers by Maddox et al. (1990a, b, c; 1996). The first ver- 
sion of the catalogue was based on 185 UKSTUr] plates with 



* United Kingdom Schmidt Telescope Unit. 



centres 8 < —20° at high Galactic latitude in the southern 
Galactic cap. The survey has since been extended to include 
the equatorial region between —17.5° < 8 < 2.5° and also 
to include equatorial regions in the northern hemisphere. 
Only the southern catalogue, as plotted in Figure 2, is used 
in this paper. Detailed analyses of the plate matching al- 
gorithm, plate matching errors, completeness, star-galaxy 
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separation and other possible sources of systematic errors 
are presented by Maddox et al. (1990 b,c; 1996). The survey 
is largely complete to bj — 20.5, though there are detectable 
systematic errors (of low amplitude) in the faint magnitude 
slice 20 < bj < 20.5. 

SFD have used the COBE/DIRBE and destriped IRAS 
maps to derive a map of the dust column density and hence 
of Galactic extinction. The Johnson B and V passbands are 
related to the APM bj passband by 

bj « B - 0.28(B - V), 

(see Maddox et al. 1990c). The SFD maps of E(B-V) can 
therefore be converted into extinction in the bj passband by 
multiplying by a factor of 4.035. The extinction computed 
from the SFD maps in the region of the southern Galactic 
pole (SGP) is plotted in Figure 2. The two plots in the lower 
panels of Figure 2 show regions of the APM Survey in which 
the extinction computed from the SFD is less than 0.2 and 
0.1 magnitudes. Evidently, Galactic extinction is relatively 
uniform and less than 0.2 magnitudes over most of the area 
of the APM survey at S < —20°. Regions of extinction higher 
than 0.2 magnitude are confined mainly to the corners at the 
top right and left of the APM map. 

Figure 3 shows dust and galaxy power spectra for var- 
ious subsets of the APM area. The power spectra in these 
figures were computed from an equal area projection, as in 
Figure 2, pixelized into 64 x 64 square pixels and applyin g an 
FFT to compute P2 (K) using the estimator of equation ( |23[ ) 
below. (These FFT estimates are not optimal, but can be 
computed very quickly. A comparison of the FFT and ML 
estimators is presented in Section 3.3.) In each panel of Fig- 
ure 3 we show power spectra for the APM Survey galaxies 
within the specified area limited at bj = 20 (filled circles). 
The crosses show power spectra for galaxies within the same 
region of sky, but with an extinction corrected magnitude, 

hf rr = bj - 4.035S(S - V), 

limited to b c ° rr < 20. (Hence the maps are regenerated 
by applying an extinction correction to each galaxy). The 
open squares show the power spectrum of the SFD extinc- 
tion map, which we have converted into modulations in the 
galaxy surface density in each pixel using 

nf* = (n)10-°- 45(4 - 03SB(fl - v) * ) . (19) 

Equation ( |l9| ) uses an approximate slope for the APM num- 
ber counts (see Maddox et al, 1990d). (Note that the mean 
extinction E(B — V)% is computed by averaging the values 
over a regular grid of 16 x 16 values within each pixel, to 
reduce the effects of small scale variations in the extinction). 

Figure 3 illustrates clearly the effects of galactic extinc- 
tion. Figure 3(a) shows that if we use the full APM survey 
area, Galactic extinction dominates the power in the APM 
Survey at wavenumbers K £ 10. Correcting the APM mag- 
nitudes for Galactic extinction results in a small reduction 
of the power at wavenumbers K ^ 10, but does not depress 
the power to the levels seen in Figures 3(c) - 3(f) for the ex- 
tinction masked APM maps. There are a number of possible 
reasons for this. The conversion from E(B~ V) to extinction 
in the bj passband may be wrong. We have tested for this 
by correlating the galaxy counts in the pixelized maps with 
E(B — V). This is plotted in Figure 4. However, as noted by 
SFD, at the limiting magnitude of the APM Survey, the fluc- 



1200 - 




600 - 



0.1 0.2 0.3 0.4 0.5 

4.035*E(B-V) 

Figure 4. Mean galaxy counts in 1.79° X 1.79° pixels plotted 
as a function of E(B — V) computed from the SFD maps. The 
error bars show the dispersion of the counts. The solid line shows 
equation ha. 



tuations in the number counts caused by galaxy clustering 
introduce a large dispersion, so it is difficult to disentangle 
the effects of Galactic extinction from galaxy clustering. The 
general trend of the counts is consistent with an extinction 
correction of 4.035-E(B — V) in the bj passband, but the 
correction is not well constrained at high extinctions. There 
may be other sources of gradients in the APM counts that 
are uncorrelated or anti-correlated with Galactic extinction, 
and so are not removed by the extinction correction. One 
effect, noted by Maddox et al. (1996) is contamination by 
stars (mainly star-galaxy mergers misclassified as galaxies) 
at low Galactic latitudes. This effect increases the counts in 
regions of high extinction, partially counteracting the effects 
of obscuration. 

In any case, simply eliminating the most highly ob- 
scured parts of the APM area has a dramatic effect on the 
power spectrum. As Figures 2 and 3 show, most of the high 
extinction regions are at S > —20° (which is why the original 
APM Survey area of Maddox et al. (1990a-c) was limited at 
this declination limit). Galactic extinction within this area 
has a negligible effect on the power spectrum except possi- 
ble at wavenumber K ,$ 3. If all pixels with 5 > —20° and 
extinctions of > 0.2 magnitudes are removed (Fig 3e), then 
the power spectrum of the extinction map is negligible at all 
wavenumbers. The power spectrum of the extinction map is 
reduced still further (by about an order of magnitude) by 
removing pixels with an extinction of > 0.1 magnitudes (Fig 
3f), whereas the power spectrum of the galaxy distribution 
hardly changes from that shown in Figures 3d-3f. This is 
powerful evidence that the power spectrum of the galaxy 
distribution in the region 8 < —20° is unaffected by Galac- 
tic extinction. In the rest of this paper, we will analyse the 
S < —20° map with a 0.2 magnitude extinction mask applied 
as in Figure 3e. 
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Figure 3. FFT power spectra for various masked cuts of the APM survey and of the SFD extinction map . The filled circles show power 
spectra for the APM Survey limited at bj = 20.0. The crosses show power spectra of APM galaxies within the same area of sky but at 
an extinction corrected magnitude limit of by rr = 20. The open squares show the power spectra of the SFD extinction maps over the 
same area (converted into a galaxy surface density modulation by equation |l9[ ). The Figures are for the following areas: (a) entire APM 
survey area as plotted in Figure 2; (b) APM area with S < —5°; (c) APM area with <5 < —5° and extinction of < 0.2 mag.; (d) APM 
area with S < -20°, (e) APM area with S < -20° and extinction of < 0.2 mag.; (f) APM area with 5 < -20° and extinction of < 0.1 
mag. 
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Figure 5. The redshift distribution of galaxies in a high com- 
pleteness subset of the 2dF Galaxy Redshift Survey. This sam- 
ple is limited at an extinction corrected magnitude limit of 
(fej) corr = 19.45. The solid line shows the fit to the dN(z) dis- 
tribution given by equation (20). 

3.2 Redshift distribution 

At the time that the BE93 and BE94 papers were written, 
few redshifts had been measured for faint galaxies in the 



APM Survey. These authors used the Stromlo/APM redshift 
survey at bright magnitudes b.j < 17 (Loveday et al, 1992) 
and the small, but deep, pencil beam surveys of Broadhurst, 
Ellis and Shanks (1988) and Colless et al. (1990, 1993) at 
bj > 20 to derive an interpolation formula for the redshift 
distribution between these magnitude limits. Here we have 
used a subset of 11120 galaxies in high completeness (> 
0.85) regions in the SGP area measured as part of the 2dF 
Galaxy Redshift Survey (2dFGRS, see Colless 1999; Folkes 
et al. 1999). The 2dFGRS uses the APM Galaxy Survey 
as the source photometric catalogue and has an extinction 
corrected (based on the SFD extinction maps) magnitude 
limit of (6j) C orr = 19.45. 

The redshift distribution for this sample is plotted in 
Figure 5. We have fitted the redshift distribution by least 
squares to a form similar to that used by BE93 

= Az 2 exp j- (±.y^dz. (20a) 

The best fitting parameters are 

z c = 0.086, P = 1.55, (bj = 19.5), (20b) 
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Figure 6. ML bandpower estimates of the two- and three-dimensional power spectrum of the APM survey limited at bj = 20.0. A 
declination limit of 8 = —20° and 0.2 magnitude extinction mask was imposed. The error bars show lcr error estimates derived from the 
Fisher matrix. The solid line shows a scale invariant linear T = 0.2 CDM model normalized so that (us)g = 1-0. 



where we have used the mean extinction correction of 0.05 
magnitudes for the 2dFGRS galaxies to convert to uncor- 
rected bj magnitudes. The fit of equation (20) is shown by 
the solid line in Figure 5. The parameters are quite close 
to those used by BE93 for bj = 19.5. To extrapola te to 
fainter and brighter magnitudes we adopt equation (20a) 
with /3 = 1.55 and adjust the parameter z c so that the me- 
dian redshift, z m , varies with magnitude limit according to 



1.36z c = 0.018(6j - 17) 15 + 0.046. 



(21) 



This formula provides an excellent fit to the median redshifts 
of published redshift surveys in the magnitude range 17 < 
b.j < 21 and to the median redshift predicted from fitting 
the lumino sity function of the 2dFGRS survey. We will use 
equations (|20a| ) and ( pl| ) to evaluate the kernel g(K/k) of 
equation (j7b(). 



3.3 Maximum likelihood power spectra of the 
APM Survey 

In this Section we show results for the maximim likelihood 
power spectra for the APM survey limited at bj = 20.0 
with a declination limit of 5 = —20° and a 0.2 magnitude 
extinction mask applied. The input maps covering the area 
shown in Figure 2 were generated with 128 x 128 pixels of 
which 4142 are 'active' (i.e. correspond to unmasked regions 
of the map). In computing the kernel g(K/k) we have as- 
sumed that the evolution parameter a — 0. The parameter 
a is not known a priori and so uncertainties in its value 
will translate into a small residual uncertainty in the am- 
plitude of the recovered three-dimensional power spectrum 



P(k), though not in its shape (see BE93; Scranton and Do- 
delson 2000). We can view the ML solution with a — 
as recovering the three-dimensional power spectrum at (ap- 
proximately) the median redshift of the survey, which for 
bj = 20 is z m « 0.14 (equation |2l|). 

The results for the two- and three-dimensional power 
spectra are shown in Figure^, together with one a error bars 
computed from the Fisher matrices. In Figure 7, we compare 
these results with those for the APM Survey limited at bj = 
19.5 and with the same sky mask. The 2d power spectrum 
for the bj — 19.5 map has a slightly higher amplitude than 
the bj = 20 estimates and is displaced slightly to higher 
wavenumbers. This is expected from the scaling properties of 
the 2d power spectrum with limiting magnitude (see BE94) . 
The 3d power spectra for the bj = 19.5 and bj = 20 maps 
are plotted in Figure (7b) and are consistent with each other. 

In Figure (8a) we compare the ML 2d power spectrum 
with the 2d power spectrum computed by applying an FFT 
to a 2048 x 2048 pixel bj = 20 APM map. The FFT power 
spectrum is computed as follows. Let ra(K) be the Fourier 
transform of the observed counts in cells of solid angle 



n(K) = J2 



me 



-iK-X| 



(22a) 



and W(JC) be the Fourier transform of the survey window 
function 



-iK Xj 



(22b) 



where Wi = 1 for active pixels and u>i = for masked pixels. 
We can form an estimate of the 2d power spectrum from 
these Fourier transforms by averaging 
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K (inverse radians) k (inverse h Mpc) 



Figure 7. ML bandpower estimates of the two- and three-dimensional power spectrum of the APM Survey. The filled circles show results 
for the APM Survey limited at bj = 20 as plotted in Figure 6. The open squares show results for a magnitude limit of bj = 19.5. For 
clarity, the bj = 19.5 points have been displaced to the right by 0.05 in the log. The lines in these figures show a linear T = 0.2 CDM 
model as in Figure 6. 




K (inverse radians) k (inverse h 1 Mpc) 



Figure 8. The filled circles show the ML bandpower estimates of P2{K) and Ps{k) as plotted in Figure 6. The open circles in Figure 
(8a) show the FFT estimates of P2(K) and lcr error bars. The open circles in Figure (8b) show the 3d power spectrum recovered by 
applying the Lucy deconvolution algorithm described by BE94 to the FFT p2{K) estimates of Figure (8a). The error bars on the Lucy 
P:i(k) estimates show the inversions obtained by fitting to the tops and bottoms of the error bars of the FFT P2{K) estimates. The solid 
line in Figure (8b) shows the scale-invariant CDM model, including a non-linear correction as described in Section 4, that provides a best 
fit to the maximum-likelihood points with fc < 0.33?iMpc _1 . The dashed line shows this model, but without the non-linear correction. 
The solid and dashed lines in Figure (8a) show the 2d power spectra computed from these fits to Pz(k). 
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ft(K) = 



|fi(K) -Are*w(K)\ 2 -jvn s 

AA 2 <?1£ K ,[#(K-K')[ 2 



(23) 



over a range of wavenumbers centred on wavenumber K. 
(We will denote this averaged estimate Pz^K)). If the aver- 
aging is done over large enough bins, so that estimates of 
P2(K) in neighbouring bins are weakly correlated, and the 
underlying fluctuations are assumed to be Gaussian, then 
the variance of P%(K) is given approximately by 



Var (A W ) = (f ) £ 



(24) 



where iVr is the total number of pixels in the map, N c is 
the number of active pixels and Nk is the number of dis- 
tinct wavenumbers used to form the average p2(K). (It is 
straightforward to derive equation (^), using the approach 
of Feldman, Kaiser and Peacock 1994.) 

The FFT power spectrum points and error bars plot- 
ted in Figure (8a) are computed from equations (^3|) and 
(pi|). These estimates are not optimal for Gaussian random 
fields, unlike the ML estimates of P2(K), and do not cor- 
rectly deconvolve the survey window function W(K). Nev- 
ertheless, the FFT power spectrum agrees well with the ML 
estimate, and even the error estimates computed from equa- 
tion (^) are in reasonable agreement with those computed 
from the Fisher matrix. At wavenumbers K £ 100, the error 
bars on the FFT estimates become smaller than the points 
on the figure, whereas the error bars on the ML estimates 
blow up. This is simply a consequence of the large pixel size 
(0 C = 0.89°) of the map used for the ML estimates compared 
to the pixel size (8 C = 3.4') of the map used for the FFT 
estimates. As explained in Section 2.3, the ML method is 
not guaranteed to be optimal or unbiased at wavenumbers 
higher than K ~ 100, since the galaxy distribution begins 
to show marked deviations from Gaussianity on these scales. 
The FFT estimator will provide an unbiased estimate at 
high wavenumbers, but the errors estimated from equation 
(p4j) will generally underestimate the true errors. 

In Figure (8b) we show the ML inversion of P-z(k) (filled 
circles) plotted against the Lucy inversion (open circles) 
of the FFT P 2 (K) estimates plotted in Figure (8a). The 
Lucy inversion algorithm used here is exactly as described 
in BE94. The error bars plotted on the open circles do not 
represent la error estimates, but simply indicate the ranges 
of the inversions found by fitting to the tops and bottoms of 
the error bars of the FFT points. The Lucy inversion is con- 
sistent with the ML estimates and, of course, extend to much 
higher wavenumbers. What is clear, however, is that the er- 
rors on the ML method at wavenumbers k £ O.l/iMpc -1 are 
large and that there is no evidence for any turnover in the 
power spectrum at smaller wavenumbers. This agrees with 
the conclusions of Eisenstein and Zaldarriaga (1999). BE93, 
BE94, Maddox et al. (1996) claimed tentative evidence for 
a turnover in the power spectrum based on an analysis of 
four separate zones of the APM Survey. However, the scatter 
in the inverted power spectra from the four zones is smaller 
than the Fisher matrix errors computed from the ML inver- 
sion. The simulations on Gaussian random fields described 
in Section (2.2) show that the Fisher matrix error estimates 
are almost certainly the more reliable. 



4 CONSTRAINTS ON CDM MODELS 

In this Section, we investigate the constraints on CDM mod- 
els from the ML estimates of Ps(k). Let P T (k) be a theoret- 
ical model for the three-dimensional power spectrum spec- 
ified by a number of parameters. We find the parameters 
that minimise 



Y^F l] (P-P T {h)) (P ] -P T (k J )), 



(25) 



where i<y is the Fisher matrix and Pi the bandpower esti- 
mates determined from the ML method. 

We first investigate simple scale-invariant (scalar spec- 
tral index n s — 1) adiabatic CDM models. These models are 
characterized by a shape parameter F and amplitude (crg) 9 . 
This amplitude may differ from the amplitude of the mass 
fluctuations, ((Ts) p , depending on the relative bias between 
fluctuations in the galaxy and the mass distributions. We 
also include a non-linear correction to the shape of the power 
spectrum using the formulae in Peacock and Dodds (1996). 
The non-linear correction requires assumptions about the 
background cosmology and the amplitude of the mass fluc- 
tuations. We fix the amplitude of the mass fluctuations to 
that inferred by Eke, Cole and Frenk (1996) from the tem- 
perature distribution of X-ray clusters: 



(<7 8 )p = 0.520, 



-0.52+0.13fi„ 



= 1, 



and adopt fl m = 0.3, Oa = 0.7, consistent with the lat- 
est results from CMB anisotropy measurements combined 
with observations of distant Type la supernovae (see e.g. 
de Bernadis et al, 2000; Hannay et al, 2000; Efstathiou et 
al, 1999). The parameters (crg) p , £l m and Oa are kept fixed 
and we vary just the two parameters (<Jg) g and T to min- 
imise x 2 ■ The non-linear corrections become significant only 
for k Si 0.2/iMpc _1 , thus provided that the sums in equa- 
tion (^i[) are restricted to low wavenumbers, the fits will be 
insensitive to the non-linear model. 

Likelihood contours for the two parameter fits are 
shown in Figures (9a) and (9b). In Figure (9a), the fit is 
restricted to wavenumbers k < 0.19/iMpc _ and in Figure 
(9b) it is resticted to k < 0.33/iMpc -1 (i.e. the first five and 
six points plotted in Figure 6a respectively). The values of 
(o~a) g and T that minimise \ 2 are very similar in these two 
cases ((o~s) g ~ 0.93, V w 0.12), but the error contours are 
much smaller for k < 0.33/iMpc" 1 . The best fitting model 
for k < 0.33/iMpc _1 is plotted in Figure 8, which also illus- 
trates the size of the non-linear correction. The error con- 
tours of Figure (9a) are probably reasonable conservative 
limits on (<Jg) g and V for the APM Galaxy Survey. Although 
the constraints on the parameters are tighter in Figure (9b), 
the wavenumber range is beginning to extend into the range 
where the non-linear correction is becoming important. 

It is not primarily the accuracy of the non-linear cor- 
rection, or its dependence on cosmological parameters, that 
makes us lean toward the more conservative limits of Fig- 
ure (9a). Rather, it is the assumption that the galaxy power 
spectrum has exactly the same shape as that of the mass 
distribution which we feel is poorly justified, especially on 
scales where the mass distribution is becoming non-linear. 
For example, Benson et al. (2000), using plausible (but phys- 
ically uncertain) assumptions about galaxy formation ap- 
plied to N-body simulations of the dark matter distribution 
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Figure 9. One, two and three cr contours for the amplitude and CDM shape parameter T determined from a likelihood fit to the APM 
three dimensional power spectrum of Figure 6 for wavenumbers (a) k < 0.19/iMpc — 1 and (b) k < 0.33/iMpc -1 . The crosses show the 
values of the parameters that minimise \ 2 - 



in a A-dominated CDM model, find that the galaxy dis- 
tribution displays non-linear biasing on scales ^ 3h~ 1 Mpc. 
In fact non-linear biasing has been found in many investi- 
gations, including those based on some of the earliest nu- 
merical simulations of the CDM model (Davis et ai, 1985). 
The nature of the bias depends on the physics of galaxy 
formation and so is difficult to predict theoretically (see 
e.g. Seljak 2000). We are therefore skeptical about using 
measurements of the galaxy power spectrum together with, 
say, CMB anisotropy measurements for the precise determi- 
nation of cosmological parameters (see e.g. Eisenstein, Hu 
and Tegmark 1998, 1999; Wang, Spergel and Strauss 1999). 
As Figure 9 demonstrates, the likelihood constraints from 
galaxy clustering are extremely sensitive to the range of 
wavenumbers used in fitting the theoretical model or (al- 
most equivalently) to the range of wavenumbers over which 
galaxies are assumed to trace the mass distribution in some 
simple way (e.g. constant bias). 

The constraints in Figure 9 are in qualitative agree- 
ment with previous analyses of the APM Galaxy Survey 
(e.g. BE93, 94; Maddox et ai, 1996) which concluded that 
the large scale clustering of the APM Survey was well fitted 
by a scale-invariant CDM model with F ~ 0.2. The con- 
tours in Figure 9a are, however, considerably tighter than 
the analogous plot in Eisenstein and Zaldarriaga (1999, their 
Figure 3). Figure (9a) is almost certainly more reliable be- 
cause it is based on a self-contained ML analysis of the APM 
map, rather than using estimates of w(8) and its errors as 
an intermediate step, as in the analysis of Eisenstein and 
Zaldarriaga. These authors argue that their limits are close 
to the (approximate) theoretical lower bounds on the pa- 
rameters (o~g) g and r computed from the formula 




4tt 



K 



8P 2 (K) 



dK, 



(26) 



Figure 10. Lower bounds on on the amplitude and shape pa- 
rameter r computed from equation (26). We show 1, 2 and 3<r 
contours for a target model CDM model with (crg) g = 1 and 
T = 0.20. 



where 5P2(K) is the difference between the true 2d power 
spectrum and a model with a different value of (o~g) g and 
F. While we agree with this equation, our evaluation of the 
constraints on (as) g and T computed from this formula differ 
somewhat from those of Eisenstein and Zaldarriaga. 

Our constraints are illustrated in Figure 10 for a target 
model with (ag) g = 1 and F = 0.20. As in Eisenstein and 
Zaldarriaga, in computing SP2(K) from -Pa(fc) we have set 
Ps(k) equal to zero for k > 0.2/iMpc _1 so that wavenumbers 
above this limit make no contribution to the \ 2 i n (M) • Our 
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Figure 11. Results of three parameter fits of CDM models to the APM 3d power spectrum. Figures 11(a) and 11(c) show one, two 
and three a likelihood contours in the (ag) g — F plane after marginalizing over n s . Figures 11(b) and 11(d) show one, two and three a 
likelihood contours in the n s — Y plane after marginalizing over (o-g) s . Figures 11(a) and (b) show results for wavenumbers k < 0.19/iMpc -1 
and Figures 11(c) and 11(d) show results for k < 0.33/iMpc" 1 . The crosses indicate the values of the parameters that maximise the 
marginalized likelihoods. 



analysis is consistent, in the sense that the error contours in 
Figure 10 are tighter than those for the real survey plotted 
in Figure (9a,). It is not clear, however, why our evaluation of 
equation (Efl) differs from that of Eisenstein and Zaldarriaga. 

We have also investigated three parameter fits to CDM 
models, varying the scalar spectral index n s in addition to 
(<rg) 9 and r. Figure 11 shows likelihood contours in the 
(as) g -T and n s -T planes after marginalizing over the third 
parameter in each case assuming a uniform prior. As in 
Figure 9, we have shown results for k < 0.19/i.Mpc -1 and 
k < 0.33foMpc _1 to demonstrate the sensitivity of the re- 
sults to the wavenumber ranges used in the fits. 

Introducing n 3 as an additional parameter significantly 
weakens the constraints on (as)g and V. It is interesting that 
from Figure (lib) we can infer, reasonably conservatively, 
that a CDM model with F « 0.5 must have a significant tilt 



of n s ^ 0.8 to be compatible with the APM power spectrum 
on large scales. 



5 CONCLUSIONS 

In this paper we have tested and applied a maximum likeli- 
hood method to estimate the 2d and 3d power spectra of the 
galaxy distribution from a two-dimensional catalogue with 
a known redshift distribution. The methods are similar to 
those applied to estimate the power spectrum C(£) of the 
CMB anisotropics. However, applied to galaxy clustering, 
the ML method provides an optimal way of estimating the 
three-dimensional power spectrum and its covariance matrix 
directly from the 2d data. This provides a simple alterna- 
tive to inverting the integral equations relating the 2d power 
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spectrum, or angular correlation function, to the 3d power 
spectrum. 

We have investigated the effects of Galactic extinction 
on the APM survey using the extinction maps of Schlegel, 
Finkbeiner and Davis (1998). Galactic extinction is shown 
to have little effect on the power spectrum over the APM 
area with 5 < —20°. Eliminating the small regions below 
this declination limit where the extinction exceeds 0.2 mag- 



nitudes! depresses the power spectrum ot dust still further so 
that it has a neghble effect on the power spectrum ot galaxy 
clustering. 

Our results show that the ML power spectra are in 
good agreement with previous estimates of the 2d and 3d 
power spectra of the APM Survey (BE93, BE94; Maddox 
et ai, 1996; Eisenstein and Zaldarriaga, 1999; Dodelson and 
Gaztahaga 2000). The ML method produces reliable esti- 
mates of the covariance matrices of the 2d power spectra. 
In agreement with Eisenstein and Zaldarriaga, we conclude 
that the errors on the 3d power spectrum have been un- 
derestimated in earlier papers (BE93, BE94) and that there 
is no evidence for a peak, or turnover, in the APM galaxy 
power spectrum at k £ 0.1/iMpc" 1 . By fitting a scale in- 
variant CDM model to the 3d power spectrum at wavenum- 
bers k < 0.19/iMpc _ we find the amplitude and shape 
parameter lie within the ranges 0.78 < (as) g < 1.18 and 
0.05 < r < 0.38 at the 2a level. Including the scalar spec- 
tral index n s as a parameter significantly weakens the con- 
strains. Nevertheless, compatibility with the APM data re- 
quires that CDM models with V ^ 0.5 have spectral indices 
n s £ 0.8 at the 2a level. 

The methods described here have applications to other 
2d surveys, for example, the forthcoming Sloan Digital Sky 
Survey (SDSS, see e.g. Margon 1999). It is also possible that 
a pixel based ML method, as described here, may prove use- 
ful for the analysis of 3d galaxy redshift surveys such as the 
2dFGRS and the SDSS. The likelihood distributions derived 
here may have some applications in parameter estimation 
studies using CMB and other data (e.g. Jaffe et al., 2000; 
Lange et al., 2000). The likelihood distributions plotted in 
Figures 9 and 11 are available from the authors on request. 
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